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Abstract 

The goal of this paper is twofold. In the first part we discuss a general approach to determine 
Lyapunov exponents from ensemble- rather than time-averages. The approach passes through the 
identification of locally stable and unstable manifolds (the Lyapunov vectors), thereby revealing 
an analogy with generalized synchronization. The method is then applied to a periodically forced 
chaotic oscillator to show that the modulus of the Lyapunov exponent associated to the phase 
dynamics increases quadratically with the coupling strength and it is therefore different from zero 
already below the onset of phase-synchronization. The analytical calculations are carried out for 
a model, the generalized special flow, that we construct as a simplified version of the periodically 
forced Rossler oscillator. 

PACS numbers: 05.45.Xt, 05.45Ac 
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I. INTRODUCTION 

As soon as synchronization phenomena in chaotic systems have been discovered 
the standard tools of nonlinear dynamics have been implemented in order to clarify this phe- 
nomenon. This is particularly true for the Lyapunov exponents (LEs), because they mea- 
sure the degree of stability along different directions and are thus the natural candidates to 
quantify the degree of synchronization of different regimes. However, several subtleties have 
been immediately discovered. For instance, the negativity of the "transversal" LE is only 
a necessary condition for the stability of complete synchronization: (i) in low-dimensional 
systems, fluctuations of the finite-time LEs may render the synchronized regime unstable 
even when the "average" exponent is negative |4 : ]; (ii) in high dimensional systems, it has 
been ascertained that the propagation of finite-am 
synchronized regime, in spite of its linear stability[5 



itude perturbations can sustain an un- 
. A still open problem concerns the 



behaviour of the LEs in the context of phase synchronization and more precisely of 

the exponent quantifying the stability of the phase dynamics. In fact, it is often claimed 
that this LE is the right order-parameter to characterize the onset of phase-synchronization: 
below the transition it is conjectured to be zero, while it is strictly negative above the tran- 



sition 



17. 



10j . However, the situation is certainly less simple than initially believed because 



a negative exponent has been found also in correspondence of locking phenomena occurring 
below the onset of phase-synchronization n| • It is therefore important to clarify analytically 
the stability of the dynamics along the "phase" direction: in the absence of coupling, this 
is a marginally stable direction and it is thus natural to expect some difficulties. Here, we 
develop a method that allows concluding that the LE corresponding to the phase dynamics 
is different from zero (and possibly positive) as soon as the coupling is switched on and 
therefore even below the onset of phase-synchronization. 

One of the main problems is the lack of analytical methods for determining even pertur- 
batively the LEs. Some ideas have been put forward for the maximum exponent jl^ . ll^ , 
because almost any initial condition eventually grows with the maximum rate and no spe- 
cial care is required to tune the direction of the perturbation. However, very little is known 
for the other exponents, starting already from the second one. This is precisely what is 
needed to determine the stability of phase dynamics in the simplest system exhibiting phase 
synchronization, i.e. in a periodically forced chaotic attractor, where the first LE accounts 



2 



for the overall instability of the chaotic dynamics. Here, we attack and solve the prob- 
lem by developing a formalism to determine LEs as ensemble- rather than time-averages. 
Similar ideas have been already discussed by Ershov and Potapov Q|, although they have 
not gone much beyond the level of formal statements. In fact, their method relies on the 
determination of the growth rates of hypervolumes of increasing dimension. While this idea 
proved very effective for the development of a powerful algorithm to compute the LEsflsj]. 
its ensemble-average extension has some limitations due to the difficulty of disentangling the 
various exponents. The advantage of our approach is that we are able to associate each non- 
degenerate LE to a field of local directions, the Lyapunov vectors (LVs). Roughly speaking, 
the ith LV is determined into two steps: the first one consists in iterating forward in time 
a hypervolume of dimension i in tangent space to identify the local orientation of the most 
expanding i directions (this is also considered in fl3|): the second step consists in iterating 
backward a vector lying within such a hypervolume. As a result, a coordinate- independent 
LV can be determined: the LE is finally obtained by averaging the corresponding instan- 
taneous expansion/contraction rate over the entire phase-space, according to the invariant 
measure. 

An objective identification of LVs is particularly interesting in the study of the hydro- 
dynamic behaviour of extended systems. In the last years, mostly as a consequence of the 
pioneering work of Posch and collaborators [tgI . \y\ . it has been discovered that in mod- 
els of fluids (more in general in Hamiltonian systems) the directions corresponding to the 
smallest (in absolute value) LEs almost coincide with long-wavelength Fourier modes. This 
observation has in turn suggested that the Lyapunov analysis naturally leads to a hydro- 
dynamic description without the need of introducing a suitable coarse graining. However, 
progress has been hindered by the lack of an absolute definition of the Lyapunov "modes" , 
that have been mostly identified with the vectors arising from the implementation of the 
Gram-Schmidt orthogonalization procedure during a standard computation of the LEs. The 



only examples of a philososphy simi 
chronotopic Lyapunov approach [If 



ar to that one outlined in the present paper concern the 



and the characterization of space-time chaos 
It is also interesting to notice that the problem of identifying the LVs is itself equivalent 
to a problem of (generalized) synchronization. In fact, the Lyapunov vectors are determined 
by integrating a skew-product system composed of the original nonlinear dynamics plus the 
"forced" evolution in tangent space. As a result, the direction of the LVs varies in a possibly 
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singular way with the position in real space. However, this difficulty does not hinder the LE 
determination, which results from an average that is substantially insensitive to the presence 
of local singularities. 

An analytic investigation of the stability of phase dynamics in a generic setup is an 
extremely difficult task because of the lack of structural stability of low- dimensional chaos. 
For this reason, it is convenient to consider suitable simplified models. The simplest system 
where phase-synchronization has been investigated is the so-called special flowj^j]. This 
is basically a skew-product system, where the phase dynamics is forced by the chaotic 
amplitude dynamics. In this system, it is possible to estabilish analytically a certain number 
of results, because the phase evolution is basically unidimensional and there is no need to deal 
with the problem of identifying the direction of perturbations. In order to perform a more 
realistic analysis of phase synchronization, a suitable coupling between phase and amplitude 
dynamics has been added t0 the special flow H.H ere ,. tt ; g up a perturbative approach 
for the weakly forced Rossler oscillator, we show that the structure of the model proposed 



in Ref. 



21( is quite similar to that one expected in generic chaotic systems, whenever the 



presence of strong dissipations allows eliminating the stable directions. Furthermore, in 



hat 



21] 



order to simplify the analytic treatment of the LEs, we focus our attention on a model t 
we call the generalized special flow (GSF), very similar to that one analyzed in Ref. 
but characterized by a finite Markov partition. As a result, we find that the modulus of the 
second LE exponent increases quadratically with the coupling strength and its corresponding 
smallness justifies the claims often found in the literature that the second LE is equal to 
zero below the onset of phase-synchronization. In other words, we conclude that the LE is 
not the right order parameter to describe this transition. 

More precisely, this paper is organized as follows. In the next section we introduce a 
general approach for the determination of Lyapunov exponents through an average over the 
invariant measure. In section ILTT1 we present our case-study model, the GSF, deriving it as a 
discrete-time approximation of a periodically forced Rossler system. In sections IIVI and El 
we illustrate the perturbation expansion for the second LE, the corresponding LV and the 
invariant measure. Finally some numerical results are presented in section IVI1 along with 
conclusions. 
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II. A GENERAL APPROACH FOR THE DETERMINATION OF LYAPUNOV 
VECTORS AND LYAPUNOV EXPONENTS 



A« = ~J rfxP(x) In 



(2) 



In this section we discuss a method to determine Lyapunov exponents from suitable 
ensemble averages. It is easy to write down a formal meaningful definition, but the problem 
lies in translating it into a workable procedure. With reference to an iV-dimensional discrete- 
time system, described by the mapping rule 

x m = f d (x t ) x g n N , (i) 

one can express the ith LE (as usual, LE are supposed to be ordered from the largest to the 
smallest one) as 

' \\djM l) (^ u2 

||V« 

where -P(x) is the corresponding invariant measure, d x fd is the Jacobian of the transfor- 
mation, and the Lyapunov vector V^(x) identifies the ith most expanding direction in 
x. 

With reference to a continuous-time system, ruled by the ordinary differential equation 

x = f c (x) x G K N . (3) 

the ith LE writes as 

n f , , \dJ c V®(x)] • V«(x) 

where • denotes the scalar product. 

Unless a clear procedure to determine the LV is given, Eqs. (|2I4|) are nothing but formal 
statements. As anticipated in the introduction, V^(x) can be obtained by following a two- 
step procedure. We start with a generic set of % linearly independent vectors lying in the 
tangent space and let them evolve in time. This is the standard procedure to determine 
LEs, and it is well known that the hypervolume Y^ identified by such vectors contains for, 
large enough times, the % most expanding directions. Furthermore, with reference to the 
set of orthogonal cordinates obtained by implementing the Gram-Schmidt procedure, the 
component of a generic vector v evolves according to the following differential equation 
(for the sake of simplicity, we refer to continuous-time systems), 



v k 



y]<7kj(x)vj l<k<i (5) 



where, as shown in Ref. [l^j], <Jkj does not explicitely depend on time, but only through 
the position x in the phase space. As a result, the zth Lyapunov exponent can be formally 
expressed as the ensemble average of the local expansion rate a iti , i.e., 

A« = f dxP(x)a M (x) (6) 



By comparing with Eq. (JU), one finds the obvious equality 

_ [d x f c V®(xj\ .V«(x) 



(7) 



||V«(x)|| 2 

In Sec. IIVI where this formalism is applied to a phase-synchronization problem, we find 
that the only workable way to obtain an analytic expression for passes through the 
determination of the direction of the corresponding LV vector V^(x). 

Let us now consider the backward evolution of a generic vector V® G Y^, Its direction 
is identified by the (i — l)-dimensional vector 

u = {ux,u 2 , • • • ,Ui-x) (8) 

where Uk = v^/vi. From the definition of u and from Eq. (j3J), one easily finds that the 
backward evolution follows the equation 

i-l 

Uk = (<?i,i ~ &k,k)uk ~ ^2 Vkjitjuj - a^i 1 < k < i (9) 

j=k+i 

This is a cascade of skew-product linear stable equations (they are stable because the Lya- 
punov exponents are organized in descending order). The overall stability is basically deter- 
mined by the smallest (&k,k ~ that is obtained for k — i — 1. It is, therefore, sufficient 
to turn our attention to the last (z — 1) component of the vector V. Its equation has the 
following structure 

ii(t) = 7 w + a(t) (10) 

where 7 = Aj — Aj_i < and we have dropped the subscript i for simplicity. The value of the 
direction u is obtained by integrating this equation. By neglecting the temporal fluctuations 
of 7 (it is not difficult to include them, but this is not important for our final goal), the 
formal solution of Eq. (|10|) reads 

u(x(t)) = / e 7( *" T) a(x) dr . (11) 



This equation does not simply tell us the value of u at time t, but the value of u when the 
trajectory sits in x(t). It is in fact important to investigate the dependence of u on x. We 
proceed by determining the deviation SjU induced by a perturbation Sxj of x along the jth 
direction, 

5 jU = ( e^-^SjG^dr (12) 

J — oo 

where, assuming a smooth dependence of o on x, (see below for a further discussion of this 
point), 

5ja(r) « a x {r)5xj{r) = a^dx^e^^ . (13) 

(notice that the dynamics is flowing backward). If the Lyapunov exponent Xj is negative, 
Sja(r) decreases for r — > — oo and the integral over r in Eq. (j!2)l converges. As a result, 
5jU is proportional to 5xj, indicating that the direction of the LV is smooth along the jth 
direction. If Xj is positive, 5jCr{r) diverges, and below time t where 

Sxjity^-^ = 1 (14) 

linearization breaks down. In this case, 5a(r) for r < to is basically uncorrelated with its 
"initial value" 5j<r(t) and one can estimate SjU, by limiting the integral to the range [t ,t] 

5ju(t) = 5xj(t) [ dre {x ^ ){t - T) a x (r) (15) 

J t 

where to is given by Eq. (fT4^) . By bounding a x with constant functions and thereby per- 
forming the integral in Eq. (|T3j) . we finally obtain 

5ju(t) w 5xj(t) + dxjity^^ (16) 

. The scaling behaviour is finally obtained as the smalles number between 1 and —j/Xj. If 
we now introduce the exponent r]j to identify the scaling behaviour of the deviation of the 
LV direction when the point of reference is moved along the jth direction in phase space, 
the results are summarized in the following way 

1 for Xj < -7 

Vj = { (17) 
— 7/Aj for Xj > —7 

The former case corresponds to a smooth behavior (the derivative is finite), while the latter 
one reveals a singular behaviour that is the signature of a generalized synchronization. 
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Although most of the assumptions made to derive the above equation are quite plausible 
(even though not rigorously proved), there is one point that needs to be more carefully 
checked: the smoothness of cx(x). In the absence of a more careful analysis of this point, 
we can only claim that the above equation provides an upper bound to the true range of 
smoothness for the LV direction. 

III. FROM THE PERIODICALLY FORCED ROSSLER SYSTEM TO THE GEN- 
ERALIZED SPECIAL FLOW 

The first model where phase synchronization has been explored is the forced Rossler 
oscillator [7J. In this section we derive a discrete-time mapping describing a forced Rossler 
system in the limit of weak coupling. We obtain what we call the Generalized Special Flow 
(GSF), because it extends a mapping previously introduced to characterize the onset of 
phase synchronization |2o| . 

The starting set of ordinary differential equations is 

x = —y — z + ey cos(f2t + ip ) 

y = x + aoy — ex sxn(Qt + ijj ) (18) 
z = a± + z(x — CI2) 

where ipo fixes the phase of the forcing term at time 0. It is convenient to introduce cylindrical 
coordinates, namely u = (if, r,z), (x = r cos <p, y = r sin 0). For the future sake of clarity, let 
us denote with S c the 3-dimensional space parametrized by such coordinates. The differential 
equation (fTH|) writes as 

u = F(u)+eG(u,ffi + W (19) 



where 



1 H — sin H — - sin 20 , a r sin 2 <p — z cos <p , a% + z(r cos <p — a 2 ) 
r 2 



T 

sin 2 cos(fit + -0>o) — cos 2 0sin(fit + -0 O ) , —j= sin20cos(f2t + -0o + vr/4) , 



F = 
G = 

Note that system (JT§j) can be written in the equivalent autonomous form 

u = F(u) +eG(u,V0, V> = 0, 



(20) 
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where if) denotes the phase of the forcing term. 

We pass to a discrete-time description, by monitoring the system each time the phase <fi is 
a multiple of 2tt. In the new framework, the relevant variables are r, z, and if), all measured 
when the Poincare section is crossed. The task is to determine the transformation mapping 
the state (r,z,if>) onto (r',z',ip'). 

In order to obtain the expression of the map, it is necessary to formally integrate the 
equations of motion from one to the next section. This can be done, by expanding around 
the unperturbed solution for e = (which must nevertheless be obtained numerically). The 
task is anyhow worth, because it allows determining the structure of the resulting map, 
which turns out to be (see appendix IA"J) 



ft = $ + (T (0) ) n + A x + e {B{ cos ip + B[ sin ip) 

r' = A 2 + e {Bl cos ip + B s 2 sin ip) (21) 

z' = A 3 + e {Bl cos ip + Bl sin ip) 

where (T^) is the average period of the unperturbed Rossler oscillator and A m 's and -B m 's 
are functions of z and r. As it is shown in appendix^! they can be numerically determined by 
integrating the appropriate set of equations. Up to first order in e, the structure of the model 
is fairly general as it is obtained for a generic periodically forced oscillator represented in 
cylindrical coordinates (as long the phase of the attractor can be unambiguously identified). 
For the usual parameter values, the Rossler attractor is characterized by a strong con- 



traction along one direction 



221 ]. As a result, one can neglect the z dependence since this 



variable is basically a function of r, and thus write 

= i/j + (T {0) ) Q + Ai(r) + e {B{{r) cos^ + B[(r) sin^) 
r > = A 2 {r) + e (B c 2 (r) cosip + B s 2 (r) sin^) (22) 

where all the functions can be obtained by integrating numerically the equations of motion 
of the single Rossler oscillator. 1 



1 Strictly speaking, A and B functions in (|21|l and arc different (see appendix El . We use the same 
notations here to simplify the presentation. 
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To simplify further manipulations, we finally recast equation (}2*2*|) in the form 

^ = ^ + K + A 1 (r)+ e gi (r) cos {ip + /3i(r)) 

r' = A 2 (r) + e g 2 (r) cos (ip + (3 2 {r)) (23) 

where 

B i( r ) = 9i(r) cos fti(r) 

B?(r) = - ft (r)sinft(r) (24) 

for i = 1,2. The parameter = (T^) VL — 2-k represents the detuning between the original 
Rossler-system average frequency and the forcing frequency Q. 

The correctness of the scheme is confirmed in Fig. ^ where all the functions defining the 
model have been numerically obtained. The very fact that they all look as one-dimensional 
curves, confirms the conjecture that ^-dependence can be neg lected. 



The GSF (J23J) generalizes the model introduced in Ref . |20j , where the effect of the phase 
on the r dynamics was not included. This implies that the GSF looses the skew-product 
structure. This has important consequences on the orientation of the second Lyapunov 
vector that we determine in the next sections. Notice also that the GSF (1231) generalizes 

n 

and justifies the model invoked in Ref. j2JJ. 

In spite of the simplification introduced by removing the z variable, a rigorous treatment 
of Eq. (}2*3*j) for generic functions g and (3 is still very difficult. A first obstacle may be the lack 
of a finite Markov partition for the unperturbed system, which does not allow us expressing 
the second order correction to the LV in a closed form (see appendix [B] for details). A second 
obstacle is that the perturbation itself may and will in general destroy the Markov parti- 
tion, making the invariant measure hardly accessible to a perturbative expansion. For both 
reasons, we restrict ourselves to considering specific A, g, and (3 functions which guarantee 
the existence of a finite Markov partitions in a finite range of the coupling constant. In the 
last section we shall comment on the possibility to extend our formalism to a more general 
setup. 

For the sake of simplicity, we have decided to analyse the following model, 



r ' = f( r ) + 2ecg(r) cos(ip + a) 
xjj' = ip + K + Ar + eb cos i/j 

10 



(25) 





FIG. 1: Numerically computed functions A;, gi and (5i (i = 1,2) for the Rossler oscillator. Rossler 
parameters have been chosen as in Ref. 2(|: ao = 0.2, a% = 1 and 02 = 9. 



where 



f(r) 



2\r\ 



g{r) 



(26) 



with r G [—1, 1]. The tent-map choice for r ensures that [—1, 0] and [0, 1] are the two atoms 
of a Markov partion. Moreover, since g(r) is equal to for r = and r = ±1, this remains 
true also when the perturbation is switched on. This is a key property that is necessary to 
perform a completely analytical treatment in the following sections. 

In this two-dimensional setup, the formal expression of the ith LE (J2J writes 



2T! 



dr / dif)P(r, ■?/>) In 
1 Jo 



|J(r,^)V ( ^(r,V;)|| 
||V w (r,^)|| 2 



(27) 
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and the Jacobian is 

/ / r (r) + 2Ecg r (r) cos(ip + a) -2ecg(r) sin(ip + a)\ 
J(r,<0) = (28) 
\ A l — sbsmif> I 

where the subscript r denotes the derivative with respect to r. The computation of the 
Lyapunov exponent therefore, requires determining both the invariant measure P(r, ip) and 
the local direction of the Lyapunov vector V^. 

IV. A PERTURB ATIVE CALCULATION OF THE SECOND LYAPUNOV EXPO- 
NENT 

In this section we derive a perturbative expression for the second LE of the GSF (|25jh 
by expanding Eq. (|27j) . One of the key ingredients is the second LV, whose direction can be 
identified by writing V = (V, 1) (for the sake of clarity, from now on, we omit the superscript 
i = 2 in V and A, as we shall refer only to the second direction). Due to the skew-product 
structure of the unperturbed map (J2*HJ) . the second LV is, for e = 0, aligned along the if) 
direction (i.e. V — 0). It is therefore natural to expand V in powers of e 

V « e«i(r,V) +e 2 v 2 (r,i)) (29) 

Accordingly, the logarithm of the norm of V is 

ln||V|| 2 = ln(l+e\ 2 ) = e 2 v\ (30) 

while its forward iterate writes as (including only those terms that contribute up to second 
order in the norm), 

jy = | £ fr(r)vi - 2ceg(r) sm(ip + a) \ 
I 1 + e(Av l - b sin ip) + e 2 Av 2 ) J 
Notice that we have omitted the (r, ip) dependence of v\ and t> 2 to keep the notation compact. 
The Euclidean norm of the forward iterate is 

1 1 JV| | 2 = 1 + 2e( At>! - b sin ip) + e 2 | (Avi - 6 sin ^) 2 + 2Aw 2 + [/ r (r)v]. - 2c^(r) sin(^ + a)] 2 j 

(32) 

and its logarithm is 

In 1 1 JV| | 2 = 2e(Av 1 - b sin if)) - e 2 1 (A^ - b sin ^) 2 - 2Av 2 - [f r (r)vi — 2cg(r) sin(ip + a)] 2 | 

(33) 
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We now proceed by formally expanding the invariant measure in powers of e 

P(r, ip) « p (ip) + ^i( r > V>) + ^(r, V)- (34) 

The determination of the pj coefficients is presented in the next section, but here we an- 
ticipate that, as a consequence of the skew-product structure for e — 0, the zeroth-order 
component of the invariant measure does not depend on the phase ip. Moreover, because 
of the structure of the tent-map, po is also independent of r, i.e. po = 1/Att. The second 
Lyapunov exponent can thus be written as 

1 



A 



1 p2ir 

dr I 

-1 Jo 



d^iy-^ — I- £Vx{ r i VO) |2e(Avi(r, ip) — b simp) — e 2 (At>i(r, ip) — bsinip) 



-2Av 2 (r, ip) + [f r (r)vi(r, ip) + 2cg(r) sm(ip + a)] 2 + v 2 (r, ip) 



+ o(e< 



(35) 



As the variable ip is a phase, it is not a surprise that some simplifications can be found 
by expanding the relevant functions into Fourier components. We start writing the first 
component of the invariant measure as 



pi(r, 



(36) 



We then turn our attention to the first order component v i(r, ip) of the second LV 1)2 9 j) . Due 
to the sinuosidal character of the forcing term in the GSF (|25|l . it is easy to verify (see the 
next section) that Vi(r,ip) contains just the first Fourier component, 



v i(r, ip) = c L(r) sm(ip + a) + R(r) cos(ip + a) 



(37) 



By now, inserting Eqs. 
obtain 



into Eq. (J33j) and performing the integration over ip, we 



9i 



L(r) sin a + i?(r) cos a 



L(r) cos a — i?(r) sin a 



b^ be r l 
+&<?!- — + A— L(r) cos a -i?(r) sin a + — (3 - A 2 ) 
8 4 L J 8 



L 2 (r) + R 2 (r) 



(38) 



+ + c — g{r)L{r) 

I r 



A/ 2 

47T 



where we have further decomposed q\{r) in its real and imaginary parts 

Qi( r ) = Qii r ) + iq\(r) 
13 



(39) 



and we have denned 



1 p2-K 

dr I dipV2(r,t()), 
-l Jo 



(40) 



which accounts for the contribution arising from the second order correction to the LV. This 
expansion shows that the highest-order contribution to the second Lyapunov exponent of the 
GSF scales quadratically with the perturbation amplitude. This is indeed a general result 
that does not depend on the particular choice of the functions used to define the GSF, but 
only on the skew-product structure of the unperturbed time evolution and on the validity 
of the expansion assumed in (|34|) (we shall comment on this last issue in the next section). 

By inserting the expression for I 2 obtained in appendix [B] (see Eq. ([Blip ) in Eq. (|38|) . we 
finally obtain the perturbative expression for the second LE, 



30 



4 



dr 



-i 



6gi(r) + -(6-A 2 ) 



L 2 (r) + R 2 (r) 



(41) 



+Acg[(r) L(r) sina + R(r) cosa 



2 r , \ t / \ Ac 2 . 

+ c — g(r)L(r) H r sin 

r 4 



Ac ( - -q\(r 



A(l 



L(r) cos a — R(r) sin a 
L(r) cos K — R(r) sin K 



Accordingly, the numerical value of the second LE can be obtained by performing integrals 
which involve the four functions q[(r), q\(r), L(r), and R(r), that are determined in the 
next section. 



V. DETERMINING THE COEFFICIENTS OF THE POWER EXPANSION 

After having more or less formally expanded the expression of the second LE in powers 
of the coupling strength e in the previous section, now we show how the coefficients can 
be determined for both the invariant measure and the direction of the LV. Notice that 
the second part of the project passes through the implementation of the general ideas put 
forward in Sec. II. 
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A. The invariant measure 



We start focusing our attention on the invariant measure P(r, if)) which can be computed 
fixed point of the Frobenius- Perron equation 

P'ir' it') = P{T ^ r) + P{r+ ^ +) (42) 
y, ^ } |detJ(r-,V>-)| |detJ(r+,V + )| 1 1 

where (r~,ip~) and (r + ,ip + ) are the two preimages of (r',if>'). It is important to notice 
that our choice of the map guarantees that two solutions do exist in the whole rectangle 
[—1, 1] x [0, 2tt] in a finite range of e- values. This will be crucial for obtained exact expressions. 
As it has been shown in the previous section, we are interested in solving the above equation 
up to first order. Accordingly, we write 

P(r', iP') = p (r', iP') + ± V ft (r)e** (43) 

Z71 z — ' 

n 

where we have expanded the first order contribution as in Eq. 1)36)1 . It is also necessary to 
expand the preimages 



r 



r J + erf (44) 



where 



^ = # + ^1 (45) 



rt = ± V" (46) 

^ = ^-K T A^- (47) 

rf = =fc — - — cos(^ + a) (48) 
1 — r' 2 

= ±cA — - — cos(^ + a) - bcosipQ (49) 



At zeroth order in e, it is easy to see that the Frobenius-Perron equation (|42|) reduces to 

Po(r',if>') = \(po(4,ip^) +Po{ro,i>o)) (50) 



whose solution is everywhere constant, as anticipated in section UVl By imposing the nor- 
malization condition, one obtains 



n = s . (51) 
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By then considering that 

| det J(r, VOI -1 = - [l + ebsm ip + ecsignr(Ag(r) sin(-?/> + a) + g r (r) cos(ip + a))] (52) 



and projecting Eq. (|42j) over its first Fourier component, we finally obtain a closed equation 
for qi(r) 

-iK 



qi{r') = - 



e- iAr °" ( 9i (r - ) - 4 - 7 e ^K") + ^ a Ag(r 



(54) 



4 4- — «" 4 

The structure of this equation is very similar to a Frobenius-Perron equation for a one- 
dimensional system. The dimensionality reduction has been made possible by exploting the 
skew-product structure of the unperturbed system. Considering also the simple expression 
of the preimages of r' (they have to be determined at zeroth order), the above equation 
can be accurately solved by implementing the standard method to solve a Frobenius-Perron 
equation (the only limit being imposed by the numerical round-off). 

In Fig. El we have plotted the real and immaginary parts of q\ for three different choices 
of the parameters b and c. In all cases, one can see a very smooth dependence, which 
thus suggests the possibility to obtain accurate fully analytic expressions by expanding 
polynomially qi{r). However, being more interested in testing the overall validity of the 
perturbative approach, we do not explore this possibility. 

In fact, in order to test the general validity of the power expansion, we have numerically 
investigated three different GSFs, corresposinding to the following choices of the functions 
/ and g: i) f(r) = 1 — 2|r|, g(r) = r 2 — \r\ as considered in (|2fij) : ii) f(r) = 0.8 — 1.8|r| 
and g(r) = 1/2, for which there is no finite Markov partition; Hi) f(r) = 2(1 — 2ec)(l — |r|) 
and g(r) = 1/2, for which the finite Markov partition existing in the unperturbed case is 
destroyed as soon as the perturbation is switched on. 

In order to compare such models, we have computed the deviation of the zero Fourier 
component of the invariant measure of the map (j23|) induced by a small finite coupling e, 

^drj^ diP[P(r,ilj)\ £ -P(r,iP)\ £=0 \. (55) 

As it can be clearly seen in Fig. in the first two cases the linear term is even equal to 
zero, while relevant multiplicative logarithmic corrections are present in the third case. This 
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FIG. 2: The first order contribution q\ to the probability density for the parameter values: A = 
—0.18, a = — 7r/4, K = 0.047T. Solid and dashed lines refer to real and imaginary parts. 

"pathological" behaviour is induced by the fact that as soon as the coupling is switched on, 
an infinite series of discontinuities in the invariant measure suddenly arises in the vicinity 
of the former fixed point r = — 1. It is, however, important to notice that no peculiarity is 
found in the more generic second case. 



B. The direction of the second Lyapunov vector 



JV 



(56) 



In this subsection we derive a self-consistent equation for the second LV. We start from 
Eq. ()31|). retaining only the relevant perturbation terms 

[frij) + 2ecg r {r) cos(^ + a)]V — 2ecg{r) sin(-?/> + a) 
AV + 1 — ebsini/j 

By computing the ratio between the components of JV we obtain the new value of the slope 

V = ev[ + ev' 2 + 

f r (r) + 2ecg r {r) cos(?/> + a) (t>i + sv 2 ) — 2cg{r) sin(-?/> + a) 



i i / 

v 1 + ev 2 



(57) 



Aevi + (1 — ebsini/j) 

where we have again kept only the relevant terms (up to first order after dividing both sides 
by e) and where v [ and v' 2 are both functions of the iterates r' and ip', 



r' = Tq + 2ecg{r) cos(^ + a) 
ip' = ip' + eb cos ij) 



(58) 
(59) 
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FIG. 3: The deviation of the zero Fourier component of the invariant measure (see Ea. (|55|) of 
map (|25j) as a function of e. Three different choices of / and g have been tested: i) f(r) = 
1 — 2\r\ and g(r) = r 2 — \r\ (crosses), ii) /(r) = 0.8 — 1.8|r| and g(r) = 1/2 (squares), and iii) 
/(r) = 2(1 — 2ec)(l — |r|) and g(r) = 1/2 (diamonds). Parameter values have been fixed to 
A = —0.18, a = — 7r/4, K = 0.04-7T, b = —0.6 and c=2.8. Abscissas are divided by e to better 
emphasize deviations from linear scaling. The logarithmic deviations dispalyed by the diamonds 
are emphasized in the inset. 

where r' = f(r) and if)' Q = ip + K + Ar are the unperturbed iterates. By replacing the 
expressions for r' and ip' in Eq. (joTj) . at leading order, we obtain 



v[(f(r),ip + K + Ar) = f r {r)vi{r,^) - 2cg(r) sin(^ + a) (60) 

As anticipated in Sec. m this recursive relation can be solved by following backwards the 
dynamics of (r,ip). It is worth stressing that the term 2cg(r) sm(tp + a), i.e. the effect of 
the phase on the amplitude dynamics, acts as a source term in Eq. (jrjU|) . In its absence, 
the latter equation would yield a trivial vanishing solution for v\ (which in turn also implies 
V2 = 0, as it can be appreciated in App. |Bj). It is therefore the feedback of the phase 
on the amplitude dynamics that generates a nontrivial orientation of the perturbed second 
Lyapunov vector. Furthermore, the structure of the source term 2cg(r) sin(ifj + a) naturally 
suggests the Ansatz (|H7|). By inserting it in Eq. (|6Tjj) we obtain two recursive equations, 



L(r) 



sign (r ) 
R(r) = 



- ( R(f(r)) sin(K + Ar) - L(f(r)) cos(K + Ar) 



9{r) 



-sign(r) [(R(f(r)) cos(K + Ar) + L{f(r)) sin(K + Ar)] 



(61) 
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FIG. 4: The functions L(r) and R(r) computed for A = —0.18 and K = (solid curve) and 
K = — 1 (dashed line). 

This equation can be solved numerically, by considering it as a recursive relation to be iter- 
ated backward in time until the fixed point solution is eventually attained and the functions 
L and R, computed with the desired precision. In Figs. H] we can see some examples of how 
they look like. 

^From the analysis carried on in Sec. |H]and in particular from Eq. ()17j) . we see that the 
condition for a smooth behaviour of the direction V along the phase-direction is (noticing 
that here, 7 = A2 — Ai) is Ai > 2A2 that is certainly verified and this fully justifies the 
expansion in Fourier modes along such a direction. On the other hand, along the expand- 
ing direction r, the codition writes A2 < 0, that is only marginally verified. The apparent 
roughness exhibited by L(r) and R(r) can therefore be a manifestation of the expected 
non-complete smoothness. It is, however, also important to stress the role played by the 
functions we have specifically considered in the GSF. In fact, the tent map induces a dis- 
continuity in the tangent space that propagates everywhere, though significantly squeezed. 
Luckily enough, as it can be appreciated in App. |Bj such singularities are integrated out 
when determining the leading contribution to the Lyapunov exponent which is therefore 
substantially insensitive and can be computed without much harm. 



19 



X/e 




-5 
-10 
-15 
-20 
-25 
-30 



0.4 
0.2 [l" 




-0.2 



"i — ' — r 



J i L 



K 1 



0.12 -0.08 -0.04 
J i I i L 



-0.1 



-0.05 







0.05 0.1 
K 



FIG. 5: Second Lyapunov exponent for the Generalized Special Flow (|25[) as a function of detuning 
K. The dashed line represent the analytical result, while dots (with the bars indicating one standard 
deviation) results of direct numerical simulations of the GSF. Parameter values have been fixed 
to: A = —0.18, a = — vr/4, b = —0.6, c = 2.8. Abscissas have been rescaled by a factor e 2 to 
evidentiate the second order coefficient. The inset magnifies a part of the larger graph. 

VI. NUMERICAL RESULTS AND CONCLUSIONS 



In this paper we have introduced a novel approach to determine analytically the Lya- 
punov exponent and applied it to the specific case of a periodically forced chaotic oscillator 
described by a model (the generalized special flow) which is also introduced here starting 
from the specific case of the Rossler oscillator. 

Given the many technical difficulties that is necessary to overcome in order to finally 
obtain the numerical value of the quadratic coefficient, it is wise to compare the analytic 
expression with the direct numerical computation of the second LE performed for small 
enough values of e. In Fig. El the second order coefficient X/e 2 is determined from the analytic 
expression (pTTJ) and by directly simulating the GSF for e- values in the range [10~ 4 ,10™ 2 ]. 
The good agreement obtained for all K values confirm the correctness of the analytical 
calculations. The relative strong negative peak around K = is a manifestation of a 
resonance phenomenon. The LE tends to be more negative when the forcing frequency is 
close to the average frequency of the chaotic attractor. 

It is also important to stress that our results are valid for arbitrarily small e, i.e. below 
the transition to phase-synchronization (if there is any) and therefore tells us that the LE 
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FIG. 6: Second Lyapunov exponent for the Generalized Special Flow (|25|) as a function of detuning 
K. Symbols and parameters values are the same as in Fig. 03 except for 6 = 0. 



corresponding to the phase dynamics is immediately different from zero, as soon as the 
coupling is switched on. 

Another important point concerns the sign of the LE: naive considerations might suggest 
that the coupling tends to stabilize the phase dynamics and thereby to give a negative LE. 
However, the left tail in Fig. El (see also the inset) definitely shows a positive exponent. It 
is desirable to find some simple heuristic arguments to understand when and why the phase 
dynamics is stable, but this does not seem to be an easy task and is left as an open problem 
for future investigations. 

The major difference between the GSF, we analyse in this manuscript and the special 
flow introduced in [^J is the term proportional to c in the equation for r in Eq. (|25jl . Such 
a term prevents the possibility of further dimension reductions and requires setting up the 
machinery we have developed in this paper. It is therefore interesting to quantify its direct 
effect on the actual value of the LE. This can be simply done, by setting the other coupling 
term h = 0, an assignment that is basically complementary to what done in the standard 
special flow. The results, reported in Fig. El show a sort of "dispersive" behaviour for the LE 
which also tends to be positive. This suggests that the back coupling of the phase dynamics 
onto the amplitude evolution maybe responsible for an eventually positive LE. 

While Eq. (14 1|) cannot by any means capture the quantitative behavior of the original 
Rossler system, still the quadratic behaviour of the second LE seems to be a very general 
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FIG. 7: Second Lyapunov exponent for the periodically forced Rdssler system close to the resonance. 
Full circles refer to a forcing frequency 0, = 1.007 (the Rdssler natural average frequency is z^o = 
1.0158(1)), while open squares correspond to f2 = 1. Rdssler parameters have been fixed as ao = 0.2, 
a\ = 1, a2 = 9, while the integration interval is about t = 10 s . The inset shows the quadratic 
relation between A^ 2 ) and e. 

feature even though we can imagine that the lack of structural stability of generic oscillators 
may mask the overall behaviour with the presence of additional stability windows. We 
have therefore computed directly the second Lyapunov exponent for the periodically forced 
Rdssler system chosing the same set of parameters (ao = 0.2, a± = 1 and a2 = 9) considered 



in Ref. 



20]. 



When the frequency of the periodic force is close to the natural frequency of the oscillator, 
Vq = 1.0158(1) for our choice of parameters, we are able to detect the quadratic behavior 
with a good accuracy, as demonstrated in Fig. Since the coupling strengths we have 
reached are much below the onset of phase synchronization, as can bee seen in Ref. (stj, 
these numerical results confirm our theoretical conclusions that the Lyapunov exponent 
corresponding to the phase dynamics deviates from zero as soon as the coupling is switched 
on. It would be now interesting to extend the analysis carried out in this paper to two coupled 
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chaotic oscillators, perhaps by investigating suitable discrete-time models such as that one 
introduced in 2^|. However, while one can presumably learn something on the phenomenon 
of phase-synchronization, we do not expect qualitative changes for the behaviour of the 
Lyapunov spectrum. This is supported by the recently reported quadratic growth of the 
fifth LE (the first one corresponding to phase dynamics) in a system of four coupled Rossler 



oscillators 
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Altogether, the numerical and analytical results presented in this paper clearly show that 
the onset of phase-sinchronization is not signalled by the second LE (or, more generally, 
the LE associated to the phase dynamics) turning neg ative, but it is rather associated to a 
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25j that is not directly related to the 



change in the structure of the dynamical attractor 
sign of the "phase exponent" . On the other hand, the quadratic dependence on the coupling 
strength makes it difficult to numerically appreciate deviations from zero (especially because 
of the statistical fluctuations that necessarily affect numerical simulations) and explains why 
in earlier studies, the LE has been mistakenly regarded as a proper order parameter to 
characterize the transition to a phase-synchronized regime. 

Another important point concerns the sign of the second Lyapunov exponent. In fact, it 
was formerly believed that phase chaos (i.e. a positive LE) can only occur in the presence 
of a specific structure of the underlying chaotic attractor (see e.g. |26[). On the other hand, 
our analytical results show that the second LE can be positive even in a context where no 
peculiar amplitude evolution has to be invoked. However, our approach does not give any 
physical insight about the expected sign of the LE. It will be certainly useful to find under 
which conditions a chaotic phase dynamics may arise. 

Finally another major achievement of this paper is that Lyapunov exponents can be ef- 
fectively determined from ensemble averages, passing through the determination of the local 
direction of the corresponding Lyapunov vectors. From a purely numerical point of view, 
there is no conceptual difficulty to applying this method for a more detailed characteriza- 
tion of high-dimensional chaos However, in the perspective of obtaining more general 
analytical results, it is desirable to go beyond systems with finite Markov partitions. 
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APPENDIX A: FROM CONTINUOUS TO DISCRETE TIME 

In this appendix we present the detailed calculations relative to the determination of the 
Poincare mapping (|21|) for the periodically forced Rossler oscillator (fTHj) . Notice, however, 
that the methodology is quite general and is indeed applicable to a generic periodically 
perturbed system, as long as it can be written in the form (|19|) . 

Let us start by introducing some useful notations 

U(r, z, ij; t) = ($(r, z, if;; t), R(r, z, if); t), Z(r, z, if>; t)) (Al) 

denotes the phase point in S c at time t of a trajectory started in (0, r, z) at time and with 
an initial phase of the forcing term equal to if) (pay attention to the fact that the triple 
(r, z, ip) G Sd)- The crossing time with the Poincare surface is determined by imposing that 
the phase $ has increased by 2n, i.e., 

$(r,z,ip;T) = 2tt. (A2) 

As we are interested in the small coupling regime, we can expand U in powers of e and 
retain just the first order term, 

U(t) = U (0) (t)+eU (1) (t) (A3) 
In particular, from Eq. (|A2|i . we obtain 



/M>(°) 

2tt = $ (0) (T) + £ $ (1) (T) = $(°)(T (0) ) +e$ (1) (T (0) ) +£--— (T (0) )T (1) (A4) 

where we have expanded T as well, assuming that T = + eT^\ Since, = 2ir, we 
conclude that 

T (i) = — L (A5) 
Ji 

etermined by the right- 

the first component of F. 



where fi = — g f T( - is determined by the right-hand side of (fT§j) with e = 0, namely, it is 
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It can be easily seen that the zeroth and first order components satisfy the differential 
equations 



Xj(o) = F(U (0)) 

U (1) = DF(U (0) )U (1) + G(U (0) ,fit + </>) 



(A6) 
(A7) 



where DF denotes the Jacobian of the velocity field F and we have introduced an explicit 
dependence on the phase ip, as it changes in going from one to the next section. The equation 
for the first order correction can be formally solved, 



U(i) 



T (0) 



drW(T<°\ r)G(U (0) (r), Qt + V) 



(A8) 



where W(i, r) is the matrix of fundamental solutions of the equation U = DFCU^ - 1 )!!. Since 
G contains only first harmonics in ip, it can be decomposed into sine and cosine components, 



G(U (0) ,fir + V) = G c (U (0) ,fir) cos^ + G S (U (0) , fir) sin^ 



where 



G c (U (0) ,fir) 



sin 2 $ (0) cos Qt - cos 2 $ (0) sin fir 
r(o) S in2$(°) cos(fir + n/4)/y/2 




(A9) 



(A10) 



and 



/ 



G s (U (0) ,fir) 



sin 2 $ (0) sin fir - cos 2 $ (0) cos fir 



\ 



-RW sin 2$(°) sm{Qr + n/4)/V2 



Accordingly, as it follows from (|A8|) . can be decomposed in the same way 



V 



(All) 



U (1) = M c (r,z) cosip + M s (r,z) smip 



where 



M c (r, z) 



T (0) 



drW(T ( V)G c (U (0) (r),nr) 



(A12) 



(A13) 



and a similar equation holds for M s (r, z). 

The first component of Eq. (jA12J) gives $W. After substituting it into Eq. ()A5|) . we 
obtain, 

r, z) cosip + M-f(r, z) sin ^) (A14) 

Ji 
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where the subscripts indicate once more the component of the vector. Accordingly, the new 
phase ip' is 

ip' = ip + QT = ip + ftT (0) (r, z) - -r-[Ml{r : z) cosip + M[(r, z) sin^] (A15) 

h 



On the other hand, from the second and third components of Eq. (|A3|) we obtain the new 
values r' and z' by also expanding the expression of T around T^°\ 

U(T) = U (0) (T (0) ) + eU (1) (T (0) ) + £ F(U (0) )T (1) (A16) 

Straightforward but tedious calculations lead to 

xj}' = ip + (T (0) )fi + A 1 + e(B{ cos tp + B{ sin ip) (A17) 
r' = A 2 +e{B^cosip + B^smip) (A18) 
z' = A 3 + e(B^cosip + B s 3 smip) (A19) 

where (T^) is the average period of the unperturbed Rossler oscillator. The functions A4 
can be determined by integrating the unperturbed equations 

A 1 (r,z) = [T {0 \r,z) - (T (0) )] fi 

A 2 (r,z) = RPXr,z',T<®) (A20) 
A 3 (r,z) = zV>(r,z;T®) 



while the functions Bf read as 



B{{r,z) = -nM[{r,z)/h 

B c 2 (r, z) = M 2 c (r, z) - f 2 M c Jh (A21) 
B c 3 (r-z) = M^z)-hM1/h 



and similar equations hold for the sine components. 



APPENDIX B: SECOND ORDER CONTRIBUTION TO THE LYAPUNOV 
VECTOR 

In this appendix we derive a closed expression for the term I 2 , which accounts for the con- 
tribution to the LE arising from second order-corrections to the LV direction (see Eq. (|40p). 
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To this pourpose, we have to consider all terms of order e in Eq. (|57|). starting from 
v[(r',ifj') = fi( r o^o) with 



Sv[ = +bc L(r' ) cos(^' + a) — R(r' ) sin(^Q + a) 



cos ip 



+2c 2 g(r) 



L r (r' Q ) sin('0Q + a) + R r (r' ) cos(ijj' + a) cos(ip + a 



(Bl) 



where r' Q = f(r) and ip' Q = ip + K + Ar are the iterates of the unperturbed GSF as defined 
in section In the following we shall not care about the possible lack of differentiability 
along the direction r for two reasons: i) we have verified that setting up a more accurate 
procedure leads to the same results, but its presentation would be more cumbersome; ii) 
the procedure is in itself correct, because in the end we are interested in the integral that is 
insensitive to the presence of singularities. 

The recursive equation for the second order term writes as 

1 



U2(r, 



fr(r) 



(B2) 



where we have introduced the source term 
1 



s(r, i/j) 



— 2cg r (r) cos(ip + a)vi(r, %p) + 5v[ + v[(r' , ip' Q )(Avi — &sin ip) (B3) 



fr(r) L 

and v[ is given by Eq. (|5Uj) . 

Being interested in the integral I2 (see Eq. (140)1 ). we see that the integration over ip can 
be easily performed since ip' ranges over [0, 2tt]. More delicate is the integral over r because 
of the folding of r' . However, one can still solve the problem by separately integrating over 
the negative and positive values of r, i.e. the two atoms of the finite Markov partition. By 
introducing the integral over negative r-values 



/■27T 

dr / dipV2(r,ip) 
1 Jo 



(B4) 



and analogously definining 11[ , we obtain from Eq. ()B2 



where S is the integral of s(r, ip). 



S 



S 



1 p2ir 

dr / dips(r,ij)) 
1 Jo 



(B5) 



(B6) 
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We thus eventually obtain 

I 2 = / 2 " + J+ = S. 
It is now convenient to express s(r-,ip) as a function of v[ only, 



(B7) 



s(r,il>) 



1 



Mr) 



bv[ sin ip + <5f ^ + 



, Av[ — 2cg r cos('0 + a) . , 



Upon integrating over ip, we obtain, 



/r(r) 



S — TTC 



i + 2c#(r)sin(^ + a)]} (B8) 



dy[L 2 (y) + i? 2 (y)] 



-27rc^ / dr (r 2 - r) sin(Ar) L r (l - 2r) cos/C - i? r (l - 2r) sinK 



+ttc 2 A / dr (r 2 - r) cos(Ar) L(l - 2r) cosiT - - 2r) sinX 

io L 

+7tc 2 / (1 - 2r) sin(Ar) L(l - 2r) cosK - R(l - 2r) sinK 

Jo L - 1 



(B9) 



After integrating by parts the integral involving L r and R r and rescaling the dummy variable, 
we finally arrive at the desired result: 



h = vrc 2 ^ j\y[L\y) + R\y)] + 



TTC 



/: 



2 / . A(l-y) 
1 ay sin ; 



L(y) cos K — sin 



(BIO) 
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